Learning Goals

Lab Description

We will work with COVID data downloaded from the New York Times. The dataset consists of COVID-19 cases and deaths in each US state during the course of the COVID epidemic.

The objective of this lab is to explore relationships between cases, deaths, and population sizes of US states, and plot data to demonstrate this

Steps

I. Reading and processing the New York Times (NYT) state-level COVID-19 data

0. Install and load libraries

1. Read in the data

cv_states_readin <- as.data.frame(data.table::fread("https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv") )


state_pops <- as.data.frame(data.table::fread("https://raw.githubusercontent.com/COVID19Tracking/associated-data/master/us_census_data/us_census_2018_population_estimates_states.csv"))

state_pops$abb <- state_pops$state
state_pops$state <- state_pops$state_name
state_pops$state_name <- NULL

cv_states <- merge(cv_states_readin, state_pops, by="state")

2. Look at the data

  • Inspect the dimensions, head, and tail of the data
  • Inspect the structure of each variables. Are they in the correct format?
dim(cv_states)
## [1] 39114     9
head(cv_states)
##     state       date fips   cases deaths geo_id population pop_density abb
## 1 Alabama 2021-08-21    1  659750  12000      1    4887871    96.50939  AL
## 2 Alabama 2021-09-17    1  764839  13048      1    4887871    96.50939  AL
## 3 Alabama 2020-07-24    1   76005   1438      1    4887871    96.50939  AL
## 4 Alabama 2021-09-03    1  714860  12394      1    4887871    96.50939  AL
## 5 Alabama 2020-03-18    1      51      0      1    4887871    96.50939  AL
## 6 Alabama 2022-02-09    1 1254032  17452      1    4887871    96.50939  AL
tail(cv_states)
##         state       date fips cases deaths geo_id population pop_density abb
## 39109 Wyoming 2021-04-09   56 56873    701     56     577737    5.950611  WY
## 39110 Wyoming 2021-09-18   56 83958    918     56     577737    5.950611  WY
## 39111 Wyoming 2020-06-01   56   910     17     56     577737    5.950611  WY
## 39112 Wyoming 2021-04-20   56 57456    705     56     577737    5.950611  WY
## 39113 Wyoming 2020-04-11   56   343      0     56     577737    5.950611  WY
## 39114 Wyoming 2021-01-06   56 45890    464     56     577737    5.950611  WY
str(cv_states)
## 'data.frame':    39114 obs. of  9 variables:
##  $ state      : chr  "Alabama" "Alabama" "Alabama" "Alabama" ...
##  $ date       : IDate, format: "2021-08-21" "2021-09-17" ...
##  $ fips       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ cases      : int  659750 764839 76005 714860 51 1254032 159169 932250 3953 104786 ...
##  $ deaths     : int  12000 13048 1438 12394 0 17452 2558 16503 114 1882 ...
##  $ geo_id     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ population : int  4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
##  $ pop_density: num  96.5 96.5 96.5 96.5 96.5 ...
##  $ abb        : chr  "AL" "AL" "AL" "AL" ...

3. Format the data

  • Make date into a date variable
  • Make state into a factor variable
  • Order the data first by state, second by date
  • Confirm the variables are now correctly formatted
  • Inspect the range values for each variable. What is the date range? The range of cases and deaths?
cv_states$date <- as.Date(cv_states$date, format="%Y-%m-%d")


state_list <- unique(cv_states$state)
cv_states$state <- factor(cv_states$state, levels = state_list)
abb_list <- unique(cv_states$abb)
cv_states$abb <- factor(cv_states$abb, levels = abb_list)


cv_states = cv_states[order(cv_states$state, cv_states$date),]


str(cv_states)
## 'data.frame':    39114 obs. of  9 variables:
##  $ state      : Factor w/ 52 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ date       : Date, format: "2020-03-13" "2020-03-14" ...
##  $ fips       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ cases      : int  6 12 23 29 39 51 78 106 131 157 ...
##  $ deaths     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ geo_id     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ population : int  4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
##  $ pop_density: num  96.5 96.5 96.5 96.5 96.5 ...
##  $ abb        : Factor w/ 52 levels "AL","AK","AZ",..: 1 1 1 1 1 1 1 1 1 1 ...
head(cv_states)
##       state       date fips cases deaths geo_id population pop_density abb
## 571 Alabama 2020-03-13    1     6      0      1    4887871    96.50939  AL
## 380 Alabama 2020-03-14    1    12      0      1    4887871    96.50939  AL
## 431 Alabama 2020-03-15    1    23      0      1    4887871    96.50939  AL
## 33  Alabama 2020-03-16    1    29      0      1    4887871    96.50939  AL
## 168 Alabama 2020-03-17    1    39      0      1    4887871    96.50939  AL
## 5   Alabama 2020-03-18    1    51      0      1    4887871    96.50939  AL
tail(cv_states)
##         state       date fips  cases deaths geo_id population pop_density abb
## 38547 Wyoming 2022-03-18   56 155907   1769     56     577737    5.950611  WY
## 39053 Wyoming 2022-03-19   56 155907   1769     56     577737    5.950611  WY
## 38906 Wyoming 2022-03-20   56 155907   1769     56     577737    5.950611  WY
## 38959 Wyoming 2022-03-21   56 155907   1769     56     577737    5.950611  WY
## 39020 Wyoming 2022-03-22   56 155988   1783     56     577737    5.950611  WY
## 38587 Wyoming 2022-03-23   56 155988   1783     56     577737    5.950611  WY
head(cv_states)
##       state       date fips cases deaths geo_id population pop_density abb
## 571 Alabama 2020-03-13    1     6      0      1    4887871    96.50939  AL
## 380 Alabama 2020-03-14    1    12      0      1    4887871    96.50939  AL
## 431 Alabama 2020-03-15    1    23      0      1    4887871    96.50939  AL
## 33  Alabama 2020-03-16    1    29      0      1    4887871    96.50939  AL
## 168 Alabama 2020-03-17    1    39      0      1    4887871    96.50939  AL
## 5   Alabama 2020-03-18    1    51      0      1    4887871    96.50939  AL
summary(cv_states)
##            state            date                 fips           cases        
##  Washington   :  793   Min.   :2020-01-21   Min.   : 1.00   Min.   :      1  
##  Illinois     :  790   1st Qu.:2020-09-05   1st Qu.:16.00   1st Qu.:  49810  
##  California   :  789   Median :2021-03-12   Median :29.00   Median : 209010  
##  Arizona      :  788   Mean   :2021-03-12   Mean   :29.78   Mean   : 536429  
##  Massachusetts:  782   3rd Qu.:2021-09-16   3rd Qu.:44.00   3rd Qu.: 665336  
##  Wisconsin    :  778   Max.   :2022-03-23   Max.   :72.00   Max.   :9079164  
##  (Other)      :34394                                                         
##      deaths            geo_id        population        pop_density       
##  Min.   :    0.0   Min.   : 1.00   Min.   :  577737   Min.   :    1.292  
##  1st Qu.:  884.2   1st Qu.:16.00   1st Qu.: 1805832   1st Qu.:   43.659  
##  Median : 3525.5   Median :29.00   Median : 4468402   Median :  107.860  
##  Mean   : 8812.6   Mean   :29.78   Mean   : 6419723   Mean   :  422.717  
##  3rd Qu.:10688.8   3rd Qu.:44.00   3rd Qu.: 7535591   3rd Qu.:  229.511  
##  Max.   :88461.0   Max.   :72.00   Max.   :39557045   Max.   :11490.120  
##                                                       NA's   :741        
##       abb       
##  WA     :  793  
##  IL     :  790  
##  CA     :  789  
##  AZ     :  788  
##  MA     :  782  
##  WI     :  778  
##  (Other):34394
min(cv_states$date)
## [1] "2020-01-21"
max(cv_states$date)
## [1] "2022-03-23"

4. Add new_cases and new_deaths and correct outliers

  • Add variables for new cases, new_cases, and new deaths, new_deaths:
    • Hint: You can set new_cases equal to the difference between cases on date i and date i-1, starting on date i=2
  • Filter to dates after July 1, 2021
for (i in 1:length(state_list)) {
  cv_subset = subset(cv_states, state == state_list[i])
  cv_subset = cv_subset[order(cv_subset$date),]

  # add starting level for new cases and deaths
  cv_subset$new_cases = cv_subset$cases[1]
  cv_subset$new_deaths = cv_subset$deaths[1]


  for (j in 2:nrow(cv_subset)) {
    cv_subset$new_cases[j] = cv_subset$cases[j] - cv_subset$cases[j-1] 
    cv_subset$new_deaths[j] = cv_subset$deaths[j] - cv_subset$deaths[j-1]
  }
  
  cv_states$new_cases[cv_states$state==state_list[i]] = cv_subset$new_cases
  cv_states$new_deaths[cv_states$state==state_list[i]] = cv_subset$new_deaths
}

# Focus on recent dates
cv_states <- cv_states %>% dplyr::filter(date >= "2021-07-01")
  • Use ggplotly for EDA: See if there are outliers or values that don’t make sense for new_cases and new_deaths. Which states and which dates have strange values?
p1 <- ggplot(cv_states, aes(x = date, y = new_cases, color = state)) +
  geom_line()

ggplotly(p1)

There were -4678 recorded new cases on January 29, 2022 in Colorado, and -4397 recorded new cases in Pennsylvania on February 8, 2022. Neither of these values seem realistic, as we cannot have negative amounts of new cases.

p2 <- ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) +
  geom_line()

ggplotly(p2)

There are a few instances in which a negative amount of people were recorded as dying on a particular day, such as -3770 deaths in Massachusetts on March 14, 2022, or -357 deaths in California on August 11, 2021.

Additionally, there are other outliers in the data where some states reported incredibly high numbers on a particular day compared to other days, such as the single-day spikes for Missouri and Tennessee.

  • Correct outliers: Set negative values for new_cases or new_deaths to 0

  • Inspect data again interactively

cv_states$new_cases[cv_states$new_cases < 0] = 0
cv_states$new_deaths[cv_states$new_deaths < 0] = 0

p3 <- ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) +
  geom_line() + geom_point(size = 0.5, alpha = 0.5)

ggplotly(p3)

5. Add additional variables

  • Add population-normalized (by 100,000) variables for each variable type (rounded to 1 decimal place). Make sure the variables you calculate are in the correct format (numeric). You can use the following variable names:

    • per100k = cases per 100,000 population
    • newper100k= new cases per 100,000
    • deathsper100k = deaths per 100,000
    • newdeathsper100k = new deaths per 100,000

  • Add a “naive CFR” variable representing deaths / cases on each date for each state

  • Create a dataframe representing values on the most recent date, cv_states_today

cv_states$per100k =  as.numeric(format(round(cv_states$cases/(cv_states$population/100000),1),nsmall=1))
cv_states$newper100k =  as.numeric(format(round(cv_states$new_cases/(cv_states$population/100000),1),nsmall=1))
cv_states$deathsper100k =  as.numeric(format(round(cv_states$deaths/(cv_states$population/100000),1),nsmall=1))
cv_states$newdeathsper100k =  as.numeric(format(round(cv_states$new_deaths/(cv_states$population/100000),1),nsmall=1))

cv_states = cv_states %>% mutate(naive_CFR = round((deaths*100/cases),2))

max_date <- max(cv_states$date)
cv_states_today = cv_states %>% filter(date==as.Date(max_date))

II. Scatterplots

6. Explore scatterplots using plot_ly()

  • Create a scatterplot using plot_ly() representing pop_density vs. various variables (e.g. cases, per100k, deaths, deathsper100k) for each state on most recent date (cv_states_today)
    • Color points by state and size points by state population
    • Use hover to identify any outliers.
cv_states_today %>% 
  plot_ly(x = ~pop_density, y = ~cases, type = 'scatter', mode = 'markers',
          color = ~state, size = ~population, sizes = c(5, 70), 
          marker = list(sizemode = 'diameter', opacity = 0.5))

Washington D.C. is a massive outlier compared to the 50 states, so we will remove it and plot the remaining areas.

  • Remove those outliers and replot.
cv_states_today %>% 
  filter(state != "District of Columbia") %>%
  plot_ly(x = ~pop_density, y = ~cases, type = 'scatter', mode = 'markers',
          color = ~state, size = ~population, sizes = c(5, 70), 
          marker = list(sizemode = 'diameter', opacity = 0.5))
  • Choose one plot. For this plot:
  • Add hoverinfo specifying the state name, cases per 100k, and deaths per 100k, similarly to how we did this in the lecture notes
  • Add layout information to title the chart and the axes
  • Enable hovermode = "compare"
cv_states_today %>% 
  filter(state != "District of Columbia") %>%
  plot_ly(x = ~pop_density, y = ~cases, type = 'scatter', mode = 'markers',
          color = ~state, size = ~population, sizes = c(5, 70), 
          marker = list(sizemode = 'diameter', opacity = 0.5),
          hover_info = "text", 
          text = ~paste0("State: ", state,
                         "<br>Cases per 100k: ", per100k,
                         "<br>Deaths per 100k: ", deathsper100k)) %>%
  layout(title = "Population-normalized cases per 100k",
         yaxis = list(title = "Cases per 100k"),
         xaxis = list(title = "Population Density"),
         hovermode = "compare")

7. Explore scatterplot trend interactively using ggplotly()

  • For pop_density vs. newdeathsper100k create a chart with the same variables using ggplotly()
  • Explore the pattern between \(x\) and \(y\)
  • Explain what you see. Do you think pop_density correlates with newdeathsper100k?
p4 <- cv_states_today %>%
  filter(state != "District of Columbia") %>%
  ggplot(aes(x = pop_density, y = deathsper100k, colour = state, size = population)) + 
  geom_point()
  
ggplotly(p4)

Population density and new deaths per 100k seem to have some amount of correlation with one another, as the deaths per 100k people in a state seems to increase as the population density of the state increases. However, the relationship seems to be non-linear, so there are likely multiple other factors at play in the relationship between these variables.

8. Multiple line chart

  • Create a line chart of the naive_CFR for all states over time using plot_ly()
    • Use the zoom and pan tools to inspect the naive_CFR for the states that had an increase in September. How have they changed over time?
  • Create one more line chart, for Florida only, which shows new_cases and new_deaths together in one plot. Hint: use add_layer()
    • Use hoverinfo to “eyeball” the approximate peak of deaths and peak of cases. What is the time delay between the peak of cases and the peak of deaths?
cv_states %>%
  plot_ly(x = ~date, y = ~naive_CFR, color = ~state, type = "scatter", mode = "lines")
cv_states %>% 
  filter(state == "Florida") %>%
  plot_ly(x = ~date, y = ~new_cases, type = "scatter", mode = "lines") %>%
  add_lines(x = ~date, y = ~new_deaths)

In many of the states, deaths seemed to peak around November 2021, whereas the new cases seemed to peak around January 2022. This is likely just because of the different strains of COVID which were rampant at these times, as the strain which was more common in January 2022 (Omicron) was significantly less deadly than earlier strains of the virus.

9. Heatmaps

Create a heatmap to visualize new_cases for each state on each date greater than July 1st, 2021 - Start by mapping selected features in the dataframe into a matrix using the tidyr package function pivot_wider(), naming the rows and columns, as done in the lecture notes - Use plot_ly() to create a heatmap out of this matrix. Which states stand out?

cv_states_mat <- cv_states %>%
  select(state, date, new_cases) %>%
  filter(date > "2021-07-01")

cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from = state,
                                            values_from = new_cases))

cv_states_mat2 <- cv_states_mat2 %>%
  column_to_rownames('date') %>%
  as.matrix()

plot_ly(x = colnames(cv_states_mat2), y = rownames(cv_states_mat2), 
        z = ~cv_states_mat2, type = 'heatmap')

The heatmap above looks almost entirely purple, and the only extraordinarily high values for cases per day seem to occur in the states with higher populations (California, New York, Texas). This information isn’t that relevant on its own, as we’d expect states with high populations to have higher case counts.

  • Create a second heatmap in which the pattern of new_cases for each state over time becomes more clear by filtering to only look at dates every two weeks (check out the omicron wave).
filter_dates <- seq(as.Date("2021-11-01"), as.Date("2022-03-22"), by="2 weeks")

cv_states_mat <- cv_states %>% select(state, date, new_cases) %>% filter( date %in% filter_dates )
cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from = state, values_from = new_cases))
rownames(cv_states_mat2) <- cv_states_mat2$date
cv_states_mat2$date <- NULL
cv_states_mat2 <- as.matrix(cv_states_mat2)

#create heatmap

10. Map

  • Create a map to visualize the naive_CFR by state on Devember 15, 2021
  • Compare with a map visualizing the naive_CFR by state on most recent date
pick.date = "2021-12-15"

# Extract the data for each state by its abbreviation
cv_per100 <- cv_states %>% filter(date==pick.date) %>% select(state, abb, newper100k, cases, deaths) # select data
cv_per100$state_name <- cv_per100$state
cv_per100$state <- cv_per100$abb
cv_per100$abb <- NULL

# Create hover text
cv_per100$hover <- with(cv_per100, paste(state_name, '<br>', "Cases per 100k: ", newper100k, '<br>', "Cases: ", cases, '<br>', "Deaths: ", deaths))

# Set up mapping details
set_map_details <- list(
  scope = 'usa',
  projection = list(type = 'albers usa'),
  showlakes = TRUE,
  lakecolor = toRGB('steelblue')
)

# Make sure both maps are on the same color scale
shadeLimit <- 125

# Create the map
# Map for today's date

# Extract the data for each state by its abbreviation
cv_per100 <- cv_states_today %>%  select(state, abb, newper100k, cases, deaths) # select data
cv_per100$state_name <- cv_per100$state
cv_per100$state <- cv_per100$abb
cv_per100$abb <- NULL

# Create hover text
cv_per100$hover <- with(cv_per100, paste(state_name, '<br>', "Cases per 100k: ", newper100k, '<br>', "Cases: ", cases, '<br>', "Deaths: ", deaths))

# Set up mapping details
set_map_details <- list(
  scope = 'usa',
  projection = list(type = 'albers usa'),
  showlakes = TRUE,
  lakecolor = toRGB('steelblue')
)

# Create the map

Deliverables

Lab 10b questions 1-2, lab 11 questions 0-9 (only first half of q9). Upload html or pdf for both lab Rmd’s to quercus.